Load data
here::set_here()
## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scripts
suppressPackageStartupMessages({
library(slingshot)
library(tradeSeq)
library(SingleCellExperiment)
library(cowplot)
library(rgl)
library(clusterExperiment)
library(RColorBrewer)
library(aggregation)
library(ggplot2)
library(pheatmap)
library(wesanderson)
library(UpSetR)
library(gridExtra)
library(scales)
})
Import data
sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
Chronological timepoints
plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[timePoint], alpha=.6)
df <- data.frame(UMAP1=reducedDim(sds)[,1],
UMAP2=reducedDim(sds)[,2],
time=factor(timePoint),
col=brewer.pal(9,'Set1')[timePoint])
ggplot(df, aes(x=UMAP1, y=UMAP2, col=time)) +
geom_point(size=.2, alpha=.5) +
scale_color_manual(values = brewer.pal(9,'Set1')) +
theme_classic() +
ggtitle("Cells colored by sampled timepoint") +
guides(colour = guide_legend(override.aes = list(alpha = 1, size=1.5)))

# ggsave("../plots/timePointDR.png", width=9)
3D plot of trajectory
plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[timePoint], alpha=.6)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=4)
plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.6)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=4)
plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.3,
box=FALSE)
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=6)
plot3d(reducedDim(sds), aspect = 'iso', col=brewer.pal(9,'Set1')[clDatta], alpha=.3,
box=FALSE, axes = FALSE, xlab = "UMAP1", ylab="UMAP2", zlab="UMAP3")
plot3d.SlingshotDataSet(sds, add=TRUE, col=c("black"), lwd=5)
axis3d('x', at = NULL, labels = FALSE, tick = TRUE, line = 0,
pos = NULL, nticks = 5)
axis3d('y', at = NULL, labels = FALSE, tick = TRUE, line = 0,
pos = NULL, nticks = 5)
axis3d('z', at = NULL, labels = FALSE, tick = TRUE, line = 0,
pos = NULL, nticks = 5)
# rgl.postscript("../plots/tmp_trajectory.pdf","pdf")
2D plots of trajectory with cell types
# pairs(reducedDim(sds), col=brewer.pal(9,'Set1')[clDatta], pch=16, cex=1/2)
dims <- c(1,2)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6),
pch=16, cex=1/2, main="Dims 1 and 2")
for(ii in seq_along(slingCurves(sds))){
c <- slingCurves(sds)[[ii]]
lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

dims <- c(1,3)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6),
pch=16, cex=1/2, main="Dims 1 and 3")
for(ii in seq_along(slingCurves(sds))){
c <- slingCurves(sds)[[ii]]
lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

dims <- c(2,3)
plot(reducedDim(sds)[,dims], col=alpha(brewer.pal(9,'Set1')[clDatta], .6),
pch=16, cex=1/2, main="Dims 1 and 3")
for(ii in seq_along(slingCurves(sds))){
c <- slingCurves(sds)[[ii]]
lines(c$s[c$ord, dims], lwd = 3, col = "black")
}

Marker genes
library(ggplot2)
markers <- c("Krt5", "Krt14",
"Sprr1a", "Trp63", #need 2 good HBC* markers here.
"Cyp2g1", "Cyp1a2",
"Ascl1", "Neurod1",
"Omp")
dims <- c(1,3)
rd <- reducedDims(sds)[,dims]
markerPlotList <- list()
for(mm in 1:length(markers)){
y <- log1p(counts[markers[mm],])
df <- data.frame(dim1=rd[,1],
dim2=rd[,2],
y=y)
markerPlotList[[mm]] <- ggplot(df, aes(x=dim1, y=dim2, colour=y)) +
geom_point(alpha=.3, size=1/2) +
scale_color_viridis_c() +
theme_classic() +
ggtitle(markers[mm]) +
theme(legend.position = "none")
}
pMarker <- cowplot::plot_grid(plotlist = markerPlotList)
pMarker

Differential expression: most significantly increasing genes in each lineage
getUpGenes <- function(yhatDf, lineage){
yhatDf1 <- yhatDf[yhatDf$lineage == lineage,]
upGenes <- yhatDf1 %>%
group_by(gene) %>%
summarize(up = yhat[20] > yhat[1]) %>%
filter(up == TRUE)
return(upGenes)
}
asTestRes <- associationTest(sce, lineages = TRUE)
## check which genes are increasing
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ tibble 3.1.0 ✔ dplyr 1.0.5
## ✔ tidyr 1.1.3 ✔ stringr 1.4.0
## ✔ readr 1.4.0 ✔ forcats 0.5.1
## ✔ purrr 0.3.4
## ── Conflicts ────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ dplyr::collapse() masks IRanges::collapse()
## ✖ dplyr::combine() masks gridExtra::combine(), Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count() masks matrixStats::count()
## ✖ dplyr::desc() masks IRanges::desc()
## ✖ purrr::discard() masks scales::discard()
## ✖ tidyr::expand() masks S4Vectors::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks S4Vectors::first()
## ✖ dplyr::lag() masks stats::lag()
## ✖ ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename() masks S4Vectors::rename()
## ✖ dplyr::slice() masks IRanges::slice()
yhatDf <- predictSmooth(sce, gene=names(sce), nPoints=20)
upGenesNeur <- getUpGenes(yhatDf, lineage=1)
upGenesSus <- getUpGenes(yhatDf, lineage=2)
upGenesHBC <- getUpGenes(yhatDf, lineage=3)
## order according to test statistic
asTestResNeurIncreas <- asTestRes[rownames(asTestRes) %in% upGenesNeur$gene,]
asTestResSusIncreas <- asTestRes[rownames(asTestRes) %in% upGenesSus$gene,]
asTestResHBCIncreas <- asTestRes[rownames(asTestRes) %in% upGenesHBC$gene,]
ooNeur <- order(asTestResNeurIncreas$waldStat_1, decreasing=TRUE)
ooSus <- order(asTestResSusIncreas$waldStat_2, decreasing=TRUE)
ooHBC <- order(asTestResHBCIncreas$waldStat_3, decreasing=TRUE)
plotlist <- list()
lineages <- c("Neur", "Sus", "HBC")
for(ll in 1:3){
curoo <- get(paste0("oo",lineages[ll]))
curres <- get(paste0("asTestRes",lineages[ll],"Increas"))
plothlp <- list()
for(kk in 1:4){
plothlp[[kk]] <- plotSmoothers(sce, counts, gene=rownames(curres)[curoo[kk]], alpha= 1/5, lwd=3/2) +
ggtitle(rownames(curres)[curoo[kk]]) +
theme(legend.position = "none")
}
plotlist[[ll]] <- plothlp
}
pNeur <- cowplot::plot_grid(plotlist=plotlist[[1]])
pSus <- cowplot::plot_grid(plotlist=plotlist[[2]])
pHBC <- cowplot::plot_grid(plotlist=plotlist[[3]])
pNeur

pSus

pHBC

## legend
pLeg <- plotSmoothers(sce, counts, gene=rownames(curres)[curoo[kk]], alpha= 1/5, lwd=3/2) +
scale_color_viridis_d(alpha = 1/5, labels = c("Neur", "Sus", "rHBC")) +
guides(colour = guide_legend(override.aes = list(alpha = 1, size=1.5)))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
pLeg

leg <- cowplot::get_legend(pLeg)
Composite main plot
p1 <- cowplot::plot_grid(NULL, pMarker, nrow=1, ncol=2, labels=c("a", "b")) #empty space for pasting 3D trajectory
p2 <- cowplot::plot_grid(pNeur, pSus, pHBC, nrow=1, labels=c("c", "d", "e"))
p3 <- cowplot::plot_grid(p2, leg, rel_widths = c(0.9,0.1))
p4 <- cowplot::plot_grid(p1,
p3,
nrow=2)
p4

ggsave("../plots/figure1.pdf", width=12, height=10)
ggsave("../plots/figure1.png", width=12, height=10)
Shared vs unique DE genes
library(UpSetR)
asTestRes2 <- tradeSeq::associationTest(sce, l2fc = log2(1.5), lineages = TRUE)
neurGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_1, "fdr") <= 0.05)]
susGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_2, "fdr") <= 0.05)]
hbcGenes <- names(sce)[which(p.adjust(asTestRes2$pvalue_3, "fdr") <= 0.05)]
geneList <- list("Neuron" = neurGenes,
"Sus" = susGenes,
"HBC" = hbcGenes)
upset(fromList(geneList), order.by = "freq")
